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Abstract 

We present the detailed analysis of the diffusive transport of spatially inhomogeneous fluid 
mixtures and the interplay between structural and dynamical properties varying on the atomic 
scale. The present treatment is based on different areas of liquid state theory, namely kinetic and 
density functional theory and their implementation as an effective numerical method via the Lattice 
Boltzmann approach. By combining the first two methods it is possible to obtain a closed set of 
kinetic equations for the singlet phase space distribution functions of each species. The interactions 
among particles are considered within a self-consistent approximation and the resulting effective 
molecular fields are analyzed. We focus on multispecies diffusion in systems with short-range 
hard-core repulsion between particles of unequal sizes and weak attractive long-range interactions. 
As a result, the attractive part of the potential does not contribute explicitly to viscosity but 
to diffusivity and the thermodynamic properties. Finally, we obtain a practical scheme to solve 
the kinetic equations by employing a discretization procedure derived from the Lattice Boltzmann 
approach. Within this framework, we present numerical data concerning the mutual diffusion 
properties both in the case of a quiescent bulk fluid and shear flow inducing Taylor dispersion. 
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I. INTRODUCTION 



Modern applications in science, medicine and technology require a better understanding 
of the molecular mechanisms controlling the flow of liquids near solid substrates and at 
interfaces [IHl]. It is well known that structural and transport properties of highly con- 
fined fluids or near free surfaces differ from their bulk behavior due to the large surface to 
volume ratio |5j. Many phenomena occurring at molecular scales such as diffusion, mixing, 
shear thinning and lane formation involve the interplay between microscopic structural and 
transport properties, which need the investigation of the long-time flow behavior. This task 
is computationally very demanding for approaches such as Molecular Dynamics, so that 
alternative methods are desirable. Some of these alternative approaches are intermediate 
between macroscopic thermodynamic and truly microscopic methods and have the scope to 
incorporate molecular details, at the price of a limited amount of numerical effort. Among 
these, the dynamic density functional theory (DDFT) and Direct Simulation Monte Carlo 
(DSMC) are prominent numerical methods. DSMC is a direct particle simulation method 
based on kinetic theory and its basic idea is to follow the trajectories of a large number of 
statistically representative particles and stochastical collisions are modeled using scattering 
probabilities. It gives results which are accurate on scales shorter than the mean free path 
[TT| [T2] . The DDFT assumes that the evolution of the system is determined by a "ther- 
modynamic force" , which is the functional derivative of the free energy functional T with 
respect the local density P-[TU]. In DDFT the state of the solute particles at time t is de- 
scribed by the average density n{r,t) while the solvent is assimilated to a continuum whose 
interactions with the solute are modeled via a stochastic heat-bath mechanism. However, 
this approach is inappropriate to describe the hydrodynamic behavior of simple liquids and 
liquid mixtures since within the DDFT picture the momentum transport can only occur via 
diffusion, but not via convection [TBHTB] . 

On the contrary, kinetic methods extending the Boltzmann equation to the dense fluid 
regime, can in principle describe both the thermodynamic and the hydrodynamic behavior of 
simple fluids. In spite of its great historical relevance in statistical physics, the Boltzmann- 
Enskog approach has rarely enjoyed the due attention in the area of inhomogeneous fluids, 
apart from some notable exceptions pTt fTS]. The reason perhaps being that, under spatially 
inhomogeneous conditions, numerical solutions of the equation are impractical. However, 
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the situation has changed with the advent of modern lattice techniques for solving the 
Boltzmann equation, collectively named the Lattice Boltzmann Method (LBM) [TnH22]- 
The simultaneous discretization of positional and translational degrees of freedom enables 
the efficient solution of such an equation by brute force. On the other hand, the application 
of the LBM to small systems is usually considered to be outside the realm of applicability 
of kinetic methods, but routinely treated within the DDFT approach, provided that the 
considered systems are not too far from local equilibrium conditions. 

In a series of recent papers we proposed a formulation of the Boltzmann-Enskog theory 
which is thermodynamically consistent, gives satisfactory values of the transport coefficients, 
and lends itself to numerical solutions within the LBM framework [231 - I25] . The method 
proved to provide reliable results in simple geometrical set-ups and was later extended to 
multicomponent fluids and to their rich and fascinating phenomenology. In the present 
paper, we investigate further issues related to the multicomponent system with special at- 
tention to the diffusive behavior. 

Following few significant studies published on the subject [SBHSU] . but differing from ours 
in the treatment of the short range correlations, we represent the evolution of the system in 
terms of the singlet phase space distribution functions, /"(r, v, t), referring to species a. The 
governing kinetic equations and the balance equations for the density and momentum current 
of the individual species, obtained in previous work employing the multicomponent extension 
of the method of Dufty and coworkers [3T], are briefly summarized in Sec. [lljto render the 
paper self-contained. These balance equations involve different kinds of forces which are 



the subject of the analysis of Sec. Ill resulting in an identiflcation of hydrostatic, capillary. 



viscous and drag forces in terms of microscopic parameters. In Sec. IV we specialize the 



theory to a binary mixture and in Sec. IV A we turn our attention to the evolution of the local 
concentration and show how to derive microscopically the advection-diffusion equation under 



suitable assumptions. In Sec. IV B we perform an hydrodynamic analysis of the coupled set 
of balance equations in order to illustrate the response of a nearly homogeneous mixture to 
small deviations from the local equilibrium state. Finally in Sec. |V] we solve numerically 
the transport equation utilizing the extension of the Lattice Boltzmann (LBM) method 
proposed in [25], where the positions are discretized on a lattice and the velocities discretized 
over a small basis set. This strategy renders the computations efficient and numerically 
stable. The method was validated against the diffusion of a small periodic inhomogeneity 
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for several values of the bulk parameters. We have also numerically studied the coupling 
between microscopic diffusion and a non uniform velocity field, a problem known as Taylor 
dispersion [32] • A numerical comparison between the analytical work and the numerical 
solution of the model shows a satisfactory agreement with the theoretical predictions. We 
conclude this section by discussing the role of the attractive tails in determining the diffusion 



coefficient. Finally, Sec, VI contains some concluding remarks. 



II. MULTICOMPONENT TRANSPORT EQUATION 

In the present paper, we shall employ a recent method to describe the isothermal transport 
properties of a mixture [25] . The idea is to simplify the transport problem by approximating 
the interaction term in such a way that non-local correlations, giving rise to the microscopic 
structure of the fluid, are taken into account. The approximation determines a non trivial 
dependence of the transport coefficients on the density profiles. In a recent paper [2S] 
we have derived the evolution of the singlet phase-space distribution function, /°'(r,v,t), 
characterizing the state of species a, of mass m", in a M-component fluid mixture, by means 
of the following transport equation: 

-r(r, V, t) + V • Vr (r, V, t) + ■ — (r, v, t) = 

-u;[r (r, V, t) - V;!(r, v, t)] + ?^ ■ (v - u(r, t))^ (r, v, t), 

(1) 

where F" is an external body force acting on species a, T the uniform temperature of the 
system and ks the Boltzmann constant. The central quantity of eq. Q is $"(r,t), which 
bears the result of collisions between particles, and whose details will be given below. In 
addition, ip"" is the local Maxwellian equilibrium of specie a, 

r(r,v,i) = n'-(r,i)[^l"^exp( ) (2) 

and the distribution shares to the same average density and velocity as the actual dis- 
tribution 

^^(r, v,t) ^ r(r, v,t){l + ^) - -(-f ■ - -(-'^)) } 

(3) 
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Eqs. (ITJ-B contain the fields n",u",u, the average partial number density of the com- 
ponent a, its average velocity and the barycentric velocity of the mixture, respectively. The 
first two quantities are defined by: 

* 'dv\ r(r,v,t). (4) 



n"(r,t)u"(r,t) I J \ v 

One also needs to specify the partial mass density, p"(r,t) = m"n"(r,t), the global number 
density, n{r,t) = J2a^°'i^j^)' global mass density 

p(r,t) = ^p"(r,t) (5) 

a 

and the barycentric average velocity at position r: 

Eq. ([T]) is an approximate isothermal representation of the revised Enskog theory (RET) 
kinetic equation [33j where, in order to obtain a workable scheme, the non-linear collision 
operator has been replaced by the two terms featuring in the r.h.s. of the equation. It 
is a simplified representation of the multicomponent RET for hard sphere mixtures, which 
contains two features that go beyond the standard Boltzmann equation approach [3l]. The 
colliding particles are separated by a distance equal to the sum of their radii and the col- 
lision frequency is modified to take into account the excluded volume effect through the 
introduction of the pair correlation function at contact in the collision integral. Such a pair 
correlation function depends on the densities through a smoothing procedure. The first term 
in the l.h.s. of eq. ([!]) describes the fast relaxation process towards local equilibrium and 
represents in an approximate fashion the non-hydro dynamic part of the collision operator. 
It contains u, a collision frequency assumed to be the same for all species. 

The form of the first term in the r.h.s of eq ([T]) is clearly reminiscent of the Bhatnagar- 
Gross-Krook (BGK) relaxation term employed in the treatment of one-component systems 
|35j . It contains an additional factor making the difference between and ip"". The factor 
multiplying the Maxwellian in eq. (|3| serves to "orthogonalize" the term —ulf"' — ip'^] to the 
term containing the effective fields, as specified below. Such a modification is necessary 
in order to produce the correct balance equation for the partial momentum and to obtain 
the correct form of the momentum equation for the individual components (see eq. (Isl) 



In the following, we consider the evolution of the partial density and of partial momentum 
current. The first is obtained by integrating eq. ([T| w.r.t. the velocity 



Ft 



p"(r,t) + V- (p°(r,t)u(r,t))+V- (p"(r,t)K(r,t)-u(r,t))=0, (7) 



where the last term in eq. ^ is the so-called dissipative diffusion current, measuring the 
drift of the a-component with respect to the center of mass velocity. 

Multiplication of eq. ([T]) by m°v and integration w.r.t. v" yields the balance equation 
for the momentum of the species a: 

^[p°(r, t)ujiv, t)] + V. (p"(r, t)<(r, t)wj(r, t) - p"(r, t)«(r, t) - w,(r, t))(wj(r, t) - w,(r, t)) 

Ffiv) $?(r,t) 
-V.vr:5.(r,t) + ^p"(r,t) + ^^p"(r,t), 

where 

7T^^ir,t) = m'' J dw{v,-Ui){v,-u,)r{v,w,t) (9) 

represents the kinetic contribution of component a to the pressure tensor. Here and in the 
following the Einstein convention on repeated indices is employed. 



III. FORCE ANALYSIS 

In ref. [23] we derived an explicit expression for the effective fields, $"(r,t), for a model 
with repulsive hard sphere potentials of different diameters, aaa and masses m", plus long 
range attractive interactions with associated potential term f/"^. The central notion is 
that this quantity is a functional of the density and velocity of each species. By treating the 
repulsive contribution in the framework of the revised Enskog theory [33j, and the attractive 
term within the random phase approximation (RPA) [36j, the effective field is represented 
as a sum of forces of different nature: 

$"(r, t) = F"'™^(r, t) + F"''^™»(r, t) + F"'^'''{r, t). (10) 

The first term represents the force acting on species a at position r due to the influence of 
all remaining particles in the system, and is the gradient of the so-called potential of mean 
force. When the system is in thermodynamic equilibrium such a force is related to the excess 
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of the chemical potential [371 EH] over its ideal gas value, fi"^^, of the a component by the 
relation: 

F"'-/(r,t) = -V/i:,,(r,t). (11) 

Explicitly, using the form of the RET collision term and an attractive potential tail, we 
obtain the following representation 

F-.-/(r, t) = -ksT J2 ^Ip I dssg^p{r, r + a^ps, t)np{r + a^ps, t) + G"'^(r, t) (12) 

where cJq^ = (cXaa + o-pp)/'^ and the integration in the first term of the r.h.s. is over the 
unit spherical surface, while the last term represents the molecular fields associated with the 
attractive forces: 

G"'^(r,t) = - j dr'n^{v\t)g^p{v,v')VrU''^{v-v') (13) 



The second and third terms of eq. (10) carry a functional dependence on the velocities, 
contributions that are neglected in semi-macroscopic models of single or multicomponents 
[39] . These terms are crucial for the correct characterization of dissipation and diffusion in 
the condensed state and result in density-dependent transport coefficients. 



The second term in the r.h.s. of eq. (10) is the drag force exerted by unlike species on 



the particle a in reason of their different drift velocities: 

F"Aa.(r,t) = - J]7"/'(r,t)(u"(r,t) -u^(r,t)), (14) 
where we have introduced an inhomogeneous friction tensor via the equation: 



7,"^(r,t) = 2(7^^ j dsSiSjgap{r,r + (T„/3S, t)n^(r + cr„/3S,t). (15) 



Finally, the last term in the r.h.s. of eq. ( 10 ) represents the viscous force due to the presence 
of velocity gradients: 



F"'--(r,t) = Y^2alJ'^-^^^^^ ! dss(7«^(r, r+a„/3§, t)^^(r+a,/3§, t)§-(u^(r+a„^s)-u^(r)), 
p V vr J 

(16) 

where gaisi'r, r + o"q,/3S, t) is the pair correlation functions at contact (|r — r'| = cxa^) and /i^/? 
is the reduced mass fiai3 = m°'+m.P colliding pair. 
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In the case of weak spatially periodic deviations from the homogeneous reference state, 
it is possible to derive explicit expressions for the effective forces discussed above. At first, 
let us consider a slowly varying periodic variation of the densities of the two species of the 
form: 

n"(r, t) = n°(t) + 5n"(t)e''i" (17) 

where qaa/s « I and are uniform densities. By substituting such a density profile into 
eq. (12) and expanding the resulting integrals up to second order in the parameter qaap, we 



find the expression: 



^ " 7 



(18) 



with w'^^ = J dr\r\"'U°'^{r). The last term in the l.h.s. of eq. (18) corresponds to the contri 



bution to the local force acting on the species a stemming from the attractive interactions 



Similarly, we estimate the viscous force by considering uniform densities and a weak 
periodic velocity field u(r,t) = u(t)e**'''", with u"^ = u'^. We find 



(r,t) ^ -|g^u,(t)e-5^a>^,„ J?^^^, (19) 



and 



F||' (r,t) ~ -— g U||(t)e'i 2^or„^<^„^W , (20) 

where we have considered the parallel and the perpendicular part of the velocity with respect 
to the wave-vector q. As a result we obtain 

= ^y/2ni^^pkBTg^pn''ali,. (21) 



IV. THE BINARY MIXTURE 

In order to proceed with analytical work it is more convenient to use as variables the local 
mass density and local momentum variables together with concentration variables. The new 
equations can be obtained by combining appropriately equations ([T]) and 

By specializing to a binary mixture, AB, the local concentration is defined as 



From the evolution equations ^ for the partial densities, the mass continuity equation reads 

9tp(r,t) + V- (p(r,t)u(r,t))=0. (23) 
and the conservation law for the local concentration 

-c(r,t) + u(r,t) ■ Vc(r,t) + -V ■ (p(r, t)c(r, t)(l - c(r,t))w(r,t))= 0, (24) 
where we have introduced the velocity difference 

w(r,t) = u^(r,t) -u^(r,t). (25) 

Using eq.([8]), the equation expressing the total momentum balance reads 

1, 
P 



dtUj{r, t) + Mi(r, i)WiUj{Y, t) + -Vivrj^^^ 



--(n^{v,t)[Ff{v) + Ff^'''f{v,t) + Ff'''''\v,t)]+n''{v,t)[Ff{^^ 0. 

(26) 



To proceed further, it is convenient to define the total local chemical potential of each species 
A{B) through the equation: 

^^P^''\-^t) - ;iA(B^V.<(^)(r,t)% - F,^(^)-^(r,t), (27) 



where we used eq. ( 11 ) for the non ideal part and the relation between the ideal gas pressure 
and the chemical potential of an ideal gas. In the isothermal system, the gradient of the 
total thermodynamic pressure is defined as 



VjP(r, t) = n^(r, t) Vj/i^(r, t) + n^(r, t) Vj/i^(r, t) (28) 

that can be seen as a special case of the Gibbs-Duhem equation. 

In the following we shall use the fact that the kinetic contribution to the gradient of the 
pressure tensor, n^^^ = ixf"- + Tif-, can be written as: 

ViTrJf ^ - 6^,V,Pu - r]^^) V.V.m, + V^m,) , (29) 

with Pid = kBT{n^{v,t) + n^{r,t)) and r/(^) = ^(n^(r,t) + n^{r,t)). As shown in Ref. 
[lO] . in the limit of small gradients also the non-ideal contribution to the viscous force in 
the momentum equation can be written as: 

J2 ^"(r, t)F"'"^^(r, t) ^ -r/(^) V^u - {-r]'^^^ + r/f V(V ■ u) (30) 

3 



where the non-ideal contribution to the shear viscosity is 
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^ ^/2TTiJ,ai3kBTall3gai3n'^n^, 



while the bulk viscosity is 



(C) _ ^„{C) 



(31) 



(32) 



Notice that, within our approximations the kinetic contribution to the bulk viscosity van- 
ishes: r]\^^ = 0. 

In order to derive an expression for w, we compute the difference between the velocities 
of the two components using eq. ([s]), and derive the following equation : 

d 



dt 



w,{r, t) + [uf (r, t)ViWj{r, t) + w,{r, t)V,uf{r, t)) 

1 rp^ip^Y \ 1 rp^ip^Y 



^f{r,t) + FHv) <l>f(r,t) + Ff(r) 



B 



(33) 



with the abbreviation u*^ = (u"^ + u-^)/2 . Before studying such an equation we shall make 
some considerations. 



A. Homogeneous diffusion 

The phenomenon of multicomponent diffusion has ever since attracted a vivid interest 
[8] . We shall specialize the discussion to the binary mixture and consider first a system 
where diffusion is the dominant mechanism to restore equilibrium, assuming that the global 
velocity of the fluid is nearly uniform. 

In eq. ([s]) the inertial term, Vj[p"^°M"], is small with respect to the terms associated 
with the viscous component of the kinetic and potential parts of the pressure tensor (see eq. 



(16) ), (ViVTj" + n^F"'^*"*"^), and their ratio is given by: 

V,[p"<M°] pu^/L puL 



n, (34) 



where L is the typical spatial scale of the gradients, u the velocity of the flow and t] the 



shear viscosity and TZ is the Reynolds number. In eq. (33) the viscous term is also negligible 
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with respect to the chemical potential term : 



(35) 



nV/i pc^L c^s ^ 

where Cs is the sound velocity and Ma = u/cs the Mach number. Thus, in the regime of low 



velocities the ratio (35) is small and the viscous force in eq. (33) can be safely neglected 



Using eq. (14), we rewrite eq. (33) as 
d 



dt 



Wi{r,t) + Vi/iD(r,i) + 



iz7r(r,t) + ^7r(r,^))-.(r,t) 



m 



B 



(36) 



where the appropriate thermodynamic field, fic, conjugated to the concentration variable, 
c = p^/ p, is the difference in the chemical potentials per unit mass of the two components 
[36] defined as: 

V,pDir,t) ^ -Lv,/x^(r,t) - -Lv,/i^(r,t). (37) 
In the homogeneous case the friction tensor is isotropic and diagonal and can be written 



as 



7 



A I a 



m 



8 y/27TpABkBT 



bulk 2 



m 



A^B 3aB ^AB 



It can also be assumed that in eq. (36) the variation in time of w is slow, so that 



w(r,t) = -^|V/i,,(r,t)- ( 



F^(r) 



F^( 



m 



B 



(38) 



(39) 



Using the Gibbs-Duhem equation, eq.(28), the chemical potential difference can be expressed 
as 

VpD 



P 



(v(/i^ - /i^) + (m^ - m^)-VP 

p 



(40) 



by substituting into eq. (36), we find that in stationary conditions, 

j^AB 



w(r,t) 



V(/i^ 



p^) + (m^ 



m^)-VP 
P 



n 



A R 

-m m 



r F^ r 



(41) 



where we have introduced the mutual diffusion coefficient D"^^ (see ref. |13]) through the 
linear relation between the velocity and the difference between the chemical potential gradi- 
ents, the factor (fc^T)"^ having been introduced in the definition for dimensionality reasons. 



By comparing eqs. (39) and (41) we find 

.AB ^bT p 



7 n m'^m^ ' 



(42) 
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which in the case of equal masses takes the simpler form D"^^ = Eq- (43) relates 



a response quantity, the friction coefficient 7 to a fluctuation quantity, D^^, the mutual 
diffusion coefficient, according to the Einstein relation. 



Relation (41) expresses the fact that the diffusion velocity is opposed to the gradient 



of the concentration field (proportional to the first term within the parenthesis) and that 



heavier molecules tend to move towards regions of higher pressure. The last term in eq. (41 ) 
corresponds to the so-called forced diffusion. We have neglected the Soret effect, that is, 
the coupling with the temperature gradient, being consistent with our isothermal treatment. 
The appropriate extension of the present theory to thermal systems was proposed in ref. 



0]. Using eq. (38), for 7, the mutual diffusion coefficient can be written explicitly as: 

D^" = A MUl (43) 



an expression identical to that derived from the Enskog analysis [H]. Moreover, assuming 
that the mass density variations are negligible and using the relation between the chemical 
potential difference hd and the Gibbs free energy per unit mass G{P,T,n,c) |19] , 

dG. 

we can write the advection-diffusion equation for the mass concentration in the suggestive 
form 

^c(r, t) + u ■ Vc(r, t) = ^ V [(c(r, t) (1 - c(r, m^^] (44) 

which bears a close resemblance with the typical DDFT equation, with the Gibbs potential 
per unit mass G replacing the Helmholtz free energy per unit volume. In the case of a binary 
ideal gas mixture {g'^^^ = 1) with equal masses, we recover the standard advection-diffusion 
equation with a constant diffusion coefficient: 

|c(r, t) + u ■ Vc(r, t) = D^^V'c{v, t) (45) 

Notice that the present diffusion coefficient D"^^ corresponds to the Enskog and not to the 
Stokes-Einstein expression, since the underlying dynamics is purely markovian [551 HU] . 

Before concluding this section, we recall that, within the random phase approximation, 
the presence of attractive tails in the pair potentials does not produce any change on the co- 
efficients of viscosity and thermal conductivity with respect to their values in the hard-sphere 
system. This result is an artifact of the RPA method and is worse than the corresponding 
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result obtained via the Enskog method [H]. However, the value of the diffusion coefficient 
does depend on the potential tails as pointed out in ref. [SO]. In fact, the diffusion current 



is proportional to the gradient of the chemical potential difference of eq. (37). 

In Fig. [l] we display the behavior of the mutual diffusion coefficient for a mixture of 
equisized hard-spheres as a function of the bulk concentration for three different values of 
the bulk packing fraction. The relative strength of the attractive tails was fixed empirically 
according to the geometric mean Lorentz-Berthelot mixing rule [51] : 



wab = V^aaWbb (46) 

and set waa = ^ksT and wbb = waa/2. 

We observe that the presence of attractive interactions tends to reduce the value of D^^ . 
The largest deviation from the unperturbed value occurs at concentration c = 1/2 and 
decreases at fixed packing fraction as the diameter increases. 



B. Hydrodynamic analysis 

We now turn our attention to the case where interspecies diffusion is coupled to acoustic 
and shear modes. The following treatment will be based on linearized equations and has 
the purpose of connecting the macroscopic hydrodynamic properties, such as the dispersion 
relations of the propagating modes and their damping, to the microscopic parameters of 



the underlying model. After linearizing eqs. (23), (24), (26) and (33) around the state 



(po, Co, u = 0, w = 0) we find the following set of equations 



dtSp{r,t)+ poV ■u{r,t) =0 

dtu{r, t) + — VP(r, t) - — (r]VMr, t) + ^7] + r]b)V{V ■ u(r, t)) 
d 

— w(r, t) + VpD(r, t) + 7w(r, t) = 
dtSc{r, t) + co(l - Co) V ■ w(r, t) = 0. 



(47) 

(48) 

(49) 
(50) 
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We now insert the trial solutions, with 5p, 5co, Uq, Wq constants, 



5p{r,t) = 5poe«*+''^"- (51) 
u{r,t) = Uoe^*+^'>-"" (52) 
6c{r,t) = fcoe^*+*'i"" (53) 
w{r,t) = woe^^+^'i-^ (54) 

(55) 



and separate the components of the velocities u and w into their longitudinal and transverse 
parts (V X u = and V • u = 0, respectively, and similarly for w.) Choosing q along the 



z-axis, we rewrite 

C5po + iqpml = (56) 

1 dP 1 dP 1/4 \ 

pQ^dp'" po^dc^p po\3 } 

(^ufy) j^q^Hufy) (58) 
Po 

C^o + '<l{^)M + + ^< = (59) 

C<(^)+7<^^^ = (60) 

C^co + iqco{l - Cq)w'q = 0, (61) 



where the upper indexes indicate Cartesian components of the vectors. We define the kine- 
matic longitudinal viscosity Ui = {Arj / 3 + rji,) / po and the kinematic shear viscosity u = i]/pq. 
Since the model is isothermal there is no coupling to the heat modes and the transverse 
velocities are completely decoupled from the remaining variables. As a consequence, the 
two shear modes describing standard diffusion of transverse momentum, can be represented 
as 

ii^(^nr,t) =MS^^^e-"^'*+^'i"- (62) 

Similarly, the transverse component of the field w decays exponentially fast due to the 
presence of internal friction 

«;=^(f)(r,t) =<^^^e-^*+''»"-. (63) 
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The remaining four longitudinal modes are mutually coupled and one has to consider the 
roots of the determinant 













iq ( dP\ 
Po V dc )P 

iqco{l - Co) C 



For the hydrodynamic analysis, it is sufficient to compute the roots of the associated fourth 
order secular equation to order g^, so to obtain the following roots: 



Cacoustic i^Cg^ 



(64) 



with a sound velocity given by 



and where 



'dp' 



, dfiD ^ ,dP^ i,dP, 



(65) 



(66) 



The last term in eq. (66) represents the damping of sound waves by interdiffusion of the 



two species. Finally the species diffusion is associated with the eigenvalue 



Cdiffv 



-Doq' 



(67) 



with 



Do = -co(l - Co)—— 
7 oc 



(68) 



which should be compared with the r.h.s of eq.(44). 



V. NUMERICAL VALIDATION 



In this section we compare some of the theoretical predictions with the numerical results 
obtained by applying the Lattice Boltzmann numerical solution of the coupled kinetic equa- 
tions Q. The discretized form of these equation has been presented in detail in appendix 
B of ref. and will not be repeated here for the sake of brevity. 
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A. Molecular diffusion 



We first determine the diffusion coefficient of an hard-sphere mixture at various packing 
fractions and compositions and then consider the effective diffusion coefficient for a system 
subject to a special type of shear flow. In the case of small perturbations around the 
equilibrium state, it is possible to obtain an analytical estimate of the so-called Taylor 
dispersion [32] . 

Let us first consider the relaxation of an initial concentration gradient in a system with 
= m^, u = and whose global density is uniform. In the initial state the composition 
varies along the z direction as a sinusoidal wave of small amplitude, A, and given by the 
two distribution functions 

/^(r,v,t = 0) = K + Asin(g,;^))e-"^^'/(2'=-^) 

/^(r,v,t = 0) = (no- Asin(g,2))e-™"'/(2fcsT)_ ^gg^ 

The diffusion constant is computed by monitoring the decay of a particular peak of ua, 
which according to the theory, decreases exponentially with an inverse characteristic time 
I/t^Qz) = D^^q^. The extracted value of D^^ as a function of the packing fraction for 
several values of the bulk composition and different diameter ratios is reported in Fig. |2] 

We observe that the mutual diffusion coefficient, D^^, increases as the concentration of 
large spheres increases at fixed value of the packing fraction, according to the theoretical 



prediction eq.(43). On the other hand, at fixed concentration and high packing fractions, 
the diffusion constant decreases as a function of the packing fraction. 

However, in the low density region we find the unexpected result that the diffusion con- 
stant increases with the packing fraction. This regime correlates with the fact that the decay 



of the perturbation (69) does not decay diffusively, but displays an oscillatory damped be- 
havior. 

This phenomenon occurs when the wavevector of the initial fluctuation is larger than a 
critical value Qc = ^J'y/AD. This apparent deviation from the standard diffusive behavior 
is the result of probing the system at small scales where standard hydrodynamics does not 
apply. However, since the phenomenon occurs only at finite wavelength below a certain 
threshold it is not in contradiction with the hydrodynamic picture presented above. The 
diffusion equation obtained in section IV A| holds in the hydrodynamic regime when, as a 
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result of many collisions, the fluid has reached local equilibrium. In terms of wave-vector 
and frequency one requires qXmf << 1 and ut << 1 , where A^/ is the mean free path and 
T the mean collision time. At densities typical of a liquid the mean free path is of the order 
of magnitude of the molecular size, while in a very diluted gas Am/ becomes large so that 
the range of validity of hydrodynamic formulae shrinks. 

The following simple analysis shows the origin of the non monotonic decay. We flrst 
decouple the "acoustic" modes in the hydrodynamic matrix, by neglecting the derivative 
of the pressure with respect to concentration, and consider the simplifled equation ( by 
neglecting ~ 0): 

C' + 7C + Ag2 = (70) 
with A = co(l — co)(^f)p- The following decay frequencies: 



display oscillatory-damped behavior for concentration fluctuations of wave-vectors q > qc, 

I 2 

with qc = y 2^- At low density we can obtain an analytic expression for such a crossover, 
since A ~ and 



o \J lixm Am/ 

where the mean free path is Am/ = ^^^^^^^2^ • terms of the wavelength = ^n/qc the 
transition from the diffusive to the oscillatory damped behavior occurs when the Knudsen 
number, expressing the ratio of the two characteristic lengths of the problem, is 

In other words, if the wavelength of the fluctuation is of the order of the mean free path, 
collisions are not frequent enough to restore local equilibrium, which is the mechanism 
determining molecular diffusion. 

Such an oscillatory decay of diffusive modes should be contrasted with the behavior 
associated to a simple BGK coUisional kernel, which does not have such oscillations [52] - [5^ . 
In fact, in the latter the friction constant, 7, is not determined self-consistently but enters 
as a free parameter and is usually assumed to be a density independent quantity. 
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B. Taylor dispersion in a periodically modulated flow 



In this subsection, we discuss a problem where one observes the interplay of a macroscopic 
and microscopic mechanisms. This occurs, for instance, when an inhomogeneous concentra- 
tion field is subjected to a non-uniform macroscopic velocity flow. As discovered by Taylor 
such a situation determines an enhancement of the molecular diffusion in the direction of 
the flow, known as Taylor dispersion j32] . 

The theoretical calculation, sketchly reported hereafter for the sake of completeness, is 
based on a multiscale perturbation analysis. We refer to the work of [HS] for mathematical 
details. 

We consider a box of length and cross-section Ly x Lz and a fluid velocity u^iy) 
periodically modulated along the y direction: 

u(r) = (^-U cos + vr^ , 0, 0^ (74) 

The boundary conditions are such that at the extremes of the box = —U, and at the 
center of the box = U. 

In the diffusive regime the concentration obeys the three dimensional advection-diffusion 



equation (45) with u given by (74). However, the description can be contracted, using 
multiscale techniques, and instead of studying the evolution of the full concentration field 
one can focus attention on the sectionally averaged concentration, C{x,t) given by: 

1 

C{x,t) = — dyc{x,y,t) (75) 



Jo 



in the presence of a laterally averaged velocity 

Uo = ^ dyu%y). (76) 
Jo 

The theory [55] shows that the salient information about the diffusive process is given by 
the following simpler one-dimensional advection-diffusion equation: 

|c(x, t) + t/o^C(x, t) = D^ff^^Cix, t) (77) 

where the new coefficient -De// is due to the renormalization of the standard molecular 
diffusion induced by the macroscopic velocity field u^{y). Its value is given by the formula: 
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where the Peclet number, Pe = ULy/2D^^ , is the ratio between the rate of advection and 
the rate of molecular diffusion. 

So far goes the theory. The above scenario can be checked with a numerical calculation 



similar to that of sub-section V A appropriately modified in order to account for the presence 
of the field u(r). We assumed an initial concentration inhomogeneity along the x direction 
under the form of two initial density fields: 

n^{x) = riQ + A sm{kxx) 

n^{x) = Uq — Asm{kxx). 

and verified that the homogeneous state is recovered exponentially with a characteristic 
time l/r^q) = D^^j-q"^, that depends on the strength of the applied velocity field and its 



wavelength as predicted by ( 78 ) . In Figs. |3j and |4| we display the different stages of 
the evolution of the concentration field in the left columns and of the velocity field in the 
right columns, obtained from our LBM code, for two different values of the strength of the 
imposed velocity field. One can see that the concentration gradient tends to decrease as the 
time increases and so does the concentration current. In a matter of ~ 30 LBM timesteps 
the concentration gradient is barely visible and the currents have faded away. During its 
evolution, the density field initially distorts in a quasi-parabolic shape and subsequently in 
a V-shaped form, being more pronounced at high Peclet. The current displays non-trivial 
patterns alternating in time and position as time proceeds. 

The numerical results obtained from our simulation are checked against theoretical pre- 
dictions, derived under the assumption that the concentration field is assimilable to a passive 
scalar. Fig. [5] displays the relaxation time for a concentration inhomogeneity for various val- 
ues of the Peclet number, and for two values of the packing fraction. The effective diffusion 
increases quadratically as a function of the Peclet number as predicted by the theory. 



VI. CONCLUSIONS 



Using a microscopic approach based on the multicomponent Boltzmann- Enskog equation 
and a self-consistent treatment of the interactions we have studied the diffusional properties 
of a mixture of hard-spheres. In order to obtain a working scheme we have employed a 
series of hypothesis and approximations. First, we have assumed that the complex many 
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body problem can be represented by means of a modified Boltzmann-Enskog equation, 
the RET, where only configurational two-particle correlations are accounted for. Since the 
RET requires a reasonable effort only in the case of hard-sphere interactions, the attractive 
potential tails have been treated within the RPA, an approximation which fails to accu- 
rately reproduce the transport coefficients. A second important simplification adopted is 
the method of Santos et al. [31], where the slowly varying, hydrodynamic fields and the 
fast non-hydrodynamic ones are decoupled at kinetic level, and the latter are treated in a 
simplified way. Third, we only considered isothermal situations for the sake of simplicity. 
The extension to non-isothermal systems will be the subject of future work. Finally, we 
have discretized the resulting transport equations on a lattice and employed the Lattice 
Boltzmann method to obtain numerical solutions. 

The present method represents a valid alternative to popular mesoscopic techniques, such 
as the pseudo-potential-Lattice-Boltzmann (LB) method or free-energy based models (see for 
instance, ref. [12] and references therein), that retain the functional form of the equilibrium 
free energy, but sacrifice the possibility of determining the transport from the microscopic 
pair potentials through controlled approximations. In contrast, our approach leads in a 
quite natural fashion to the determination of thermodynamic forces compatible with the 
free energy methods, but in addition determines self-consistently the non equilibrium forces 
necessary to guarantee the correct hydrodynamic behavior. 

We have obtained a derivation of the advection-diffusion equation for the concentration 
and the self-consistent determination of the diffusion coefficient, which in the homogeneous 
case reduces to the Chapman-Enskog value. The study of the long wavelength and low 
frequency properties of the model has been performed and agrees with the results obtained 
by standard hydrodynamic analysis [361 EH] of mixtures. 

A second merit of the present formulation is to lend itself to numerical solution via the 
Lattice Boltzmann method. Our computational approach takes into account the dynamics 
of flowing liquids on space-time scales of hydrodynamic interest. These scales are out of 
reach for Molecular Dynamics, which in principle describes ab-initio the system, since the 
probabilistic nature of the singlet distribution function does not require averaging the data 
as in particle-based methods. In addition, the proposed method can cope very naturally 
with the critical situations of low concentrations of one species. 

By simple numerical experiments, we have verified that the present version of the LBM 
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allows to extract the value of the diffusion coefficient from the decay of small periodic con- 
centration fluctuations. Moreover, we considered cases where the Taylor dispersion mecha- 
nism provides an enhancement of diffusion, thus further showing that the present numerical 
scheme is capable of handling molecular mechanisms together with driving forces acting on 
much larger scales. 

We plan future applications of the present approach to the study of non-uniform sub- 
strates, multiphase flows and transport in narrow channels. 
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FIG. 1: Ratio between the mutual diffusion coefficient with attractive tails and the coefficient of 
a system without attractive tails as a function of concentration. The mixture consists of equisized 
spheres with attractive potentials whose strength is chosen according to the Lorentz-Berthelot 



mixing rule (46). The hard sphere radii are uaa = = 2 in one case and gaa = = 4 and 
the packing fraction ^3 = ^{n^cr\A + '>T'^(^bb) kept fixed at values 0.6 and 0.3, while varying 
concentration. The effect of the potential tails is the largest for equal concentrations. 
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FIG. 2: Numerical test of the diffusion process in bulk conditions. The vertical axis represents 
a measure of the mutual diffusion coefficient obtained from LBM simulations (all data expressed 
in LBM timestep units). We monitored the evolution towards equilibrium of a sinusoidal concen- 



tration fluctuation (see eq. (69)) of wave- vector and extracted the characteristic decay time, 
1/T{qz) = D'^^q^. The plots report the inverse decay time versus packing for various values of the 
composition and diameter ratio and for a fixed value of = 40. 
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FIG. 3: Time evolution of the density (left column) and current (right column) in Taylor dispersion. 
The initial concentration modulation is along the x direction whereas the external field varies along 
the y direction according a cosine law. In the left column we report the evolution of the density 
of the large species every 5 LBM timesteps. In the right column we report the evolution of the 
associated current in the y direction. Data correspond to aAA = 8, asB = 4, cq = 0.5, Pe = 1, 
average packing ^3 = 0.211, and for a simulation box of 80 x 40 x 40. The color scale refers to both 
the density and current plots. Both reported data are normalized according to the initial values of 
the respective fields. 
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FIG. 5: Mutual diffusion coefficients in Taylor dispersion obtained for packing fraction of 0.211 
(circles) and 0.332 (squares respectively). Filled symbols correspond to a box of length = 80, 



while open symbols to a box of length L^^ = 40. The dashed line corresponds to eq. (78). 
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